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Accurate color reproduction is important in many applications of 3D print¬ 
ing, from design prototypes to 3D color copies or portraits. Although full 
color is available via other technologies, multi-jet printers have greater po¬ 
tential for graphical 3D printing, in terms of reproducing complex appear¬ 
ance properties. However, to date these printers cannot produce full color, 
and doing so poses substantial technical challenges, from the shear amount 
of data to the translucency of the available color materials. In this paper, we 
propose an error diffusion halftoning approach to achieve full color with 
multi-jet printers, which operates on multiple isosurfaces or layers within 
the object. We propose a novel traversal algorithm for voxel surfaces, which 
allows the transfer of existing error diffusion algorithms from 2D printing. 
The resulting prints faithfully reproduce colors, color gradients and fine- 
scale details. 

Categories and Subject Descriptors: 

General Terms: 3D color printing, half-toning, error diffusion 


1. INTRODUCTION 

In this paper, we present algorithms for producing 3D color prints, 
which are highly accurate and detailed. From a technical stand¬ 
point, we consider the problem of halftoning a signal defined on a 
surface for the purposes of accurate color reproduction in 3D print¬ 
ing when using translucent materials. More specifically, we present 
algorithms for error-diffusion on voxel representations of surfaces, 
which are compatible with translucent printing materials. 

Accurate color reproduction is important in many applications of 
3D printing, especially for design-prototypes or 3D color copies, 
where a texture-mapped 3D scan of an object is to be reproduced 
in a color-consistent way; this includes 3D portraits. 

Whi le 3D color printin g has existed for som e time using powder- 
binde r | |3DSystems 2014| and layer-laminated |MCor Technologies] 
|2014) systems, multi-jet 3D printers with sufficient materials for 
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color printing have only recently come on the market |Stratasys] 
|2014| . In comparison to powder-based systems, multi-jet printers 
produce smoother surfaces at higher resolutions, and allow more 
control over the internal structure of the print. In the longer term, 
such printers have much greater potential to reproduce complex ap¬ 
pearance properties. In fact, multi-jet printers with fewer materials 
have been used to approximate desired single-tone sub surface scat¬ 
tering properties [Hasan et al. 2010[|Dong et al. 2010|, a nd opaque 
objects inside transparent shells [Vidimce et al. 2013) . Unfortu¬ 
nately, software allowing full color printing with multi-jet printers 
is not yet available; the factory software allows colors to be selected 
from a palette of 50 and manually assigned to shells. 

Without accurate color reproduction, more complex appearance 
characteristics are limited in their applicability and appeal. There¬ 
fore, accurate color reproduction is a crucial step to fully exploiting 
the appearance capabilities of these printers. We provide precisely 
this critical building block for graphical 3D printing with this paper. 

However, there are substantial technical challenges in color re¬ 
production with multi-jet printers. First, is the amount of data that 
must be processed. For accurate color, the material assigned to each 
voxel must be controlled. Modern 3D printers have a resolution of 
more than 18 million voxels per cubic centimeter. Thus, prints such 
as those in Figure [^require tens of billions of voxels (see Table |I| 
for details). Both resolution and build volume of 3D printers are 
expected to increase in the coming years. 

Second, for multi-jet printers currently on the market, the colored 
materials are highly translucent, which means that the organization 
of the materials a significant distance beneath the surface, up to a 
few millimeters, can greatly affect the perceived color as discussed 
in Section!^ This translucency leads to both blurring of fine-scale 
details in textures, and to artifacts caused by severe dot-gain if care 
is not taken in the material arrangement. It is therefore highly im¬ 
portant to control material placement several layers of voxels be¬ 
neath the surface, considering the transmission and scattering prop¬ 
erties of the materials, greatly complicating the computational as¬ 
pects of halftoning algorithms. 

While we expect less translucent materials to come on the mar¬ 
ket in the future, not only will our proposed technique also work 
for these materials, it will actually perform better-producing equal 
or larger color gamuts with less computational effort. Further, em¬ 
ploying translucent materials may have its own advantages, when 
we consider materials with similar translucency as observed in skin. 

Third, while powder-binder and layer-laminated systems, like 
2D color printing, may combine up to 3 inks at a single surface 
location, multi-jet printing allows exactly one material per voxel. 

In this paper, we leverage the knowledge of decades of re¬ 
search in color imaging, color management and 2D color printing, 
to maximize the quality and exploit the full capabilities of high- 
resolution multi-material 3D printers-and push their limits towards 
realism. To this end, we propose a geometry-adaptive error diffu¬ 
sion halftoning algorithm, which includes the following technical 
contributions: 
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Fig. 1. Color 3D prints computed with our software and printed with a multi-material printer show tremendous details and realism. Three of the apples and 
the thumb are not printed. 


—A traversal algorithm for voxel representations of surfaces, 
which maps 2D anisotropic error diffusion filters onto a surface 
in a consistently oriented way (Section|^ and requires only local 
information to do so. 

—A layered halftoning algorithm, which combines the traversal al¬ 
gorithm with an arbitrary 2D error diffusion algorithm, and can 
adapt to the translucency of the materials or increase the color 
gamut by varying the number of layers (Sectionj^. 

We do not attempt to approximate the full BSSRDF of an object, 
but rather present halftoning algorithms capable of reproducing an 
object’s albedo texture even when the printing materials are highly 
translucent, under non-pathological illuminaton conditions, such as 
a light source in contact with (as in Figure]^ or inside the printed 
object. Note that a high-quality halftone, such as those produced by 
our algorithms, is necessary, but not sufficient to accurately repro¬ 
duce an object’s color. Additionally, one needs colorimetrically or 
spectrally accurate measurements of the object’s color, which are 
typically not provided by common 3D scanners, and an accurate 
optical printer model. In Section |7.2| we describe a fully empiri¬ 
cal model for our printer. A more accurate model is left for future 
work. 

2. RELATED WORK 

Existing 3D color printing tech nologies. Powder-binder |3DSys-| 
Items 20l^ and layer-laminated |MCor Technologies 20i4| sys- 
tems provide full color solutions. In both cases, CMYK inks are 
applied to a white base material. Although these technologies are 
effective for color printing, they do not provide the resolution or 
smooth surfaces of multi-jet systems and are limited to nearly 
opaque materials. They therefore have less potential for future com¬ 
binations of color with complex appearance properties. 
Appearan ce reproduction and multi-material fabrication. 
Spec2Fab | |Chen et al. 2013) is a general reducer-tuner framework 
for specification-driven digital fabrication, which allows textures 
to be replicated on a 3D print. They use an error diffusion opti¬ 
mization of material layerings, effectively a contone, with a uni¬ 
form error filter (error is pushed equally to all neighbors). Although 
an important first step in texture-mapping for multi-material print¬ 
ers, it does not allow for anisotropic error diffusion filters, and 
the i terative optimization p rohibits a streaming architecture. Open- 
Fab |Vidimce et al. 2013) is a programmable fabrication pipeline 
for 3D printers, which uses in-slice 2D error diffusion dither¬ 
ing [Floyd and Steinberg 1976| . The authors observe that by dither¬ 
ing in 3D they could avoid streaks. Our approach treats the color 


signal where it is defined-on the surface-by mapping 2D filters into 
the tangent space of the surface. As discussed below, this allows 
us to colorimetrically characterize the 3D printer in a geometry- 
independent way. 

Multi-material 3D Printing has b een used to reproduce speci¬ 
fied su bsurface scattering properties [Hasan et al. 2010|[Dong et ah] 
|2010| . Fabrication of directional BRDFs for planar or near-planar 
surfac es has been done usi ng multi-material printing [Lan et ah] 
|2013) and photolithography [Levin et al. 2013| . 

Given the work that has already been done using multi-material 
3D printing to reproduce complex appearance properties, the intro¬ 
duction of full color opens up tremendous possibilities for future 
appearance fabrication. 

Tone reproduction in FDM prints. Some recent work has focused 
on the challenging task of improving tone reproduction in f used de- 
position modeling (FDM) printing. Hergel and Lefebvre [Hergelj 
land Lefebvre 2014| optimize seam placement in multi-filament 
FDM prints to hide or reduce artifacts from changing filaments. 
Reiner et al. [Reiner et al. 2014) perform a type of halftoning for 
FDM printers while maintain long filament paths. Switching fila¬ 
ment heads not only creates artifacts, but also increases print time. 
Both of these methods are specific to FDM printers. 

3D Halftoning. 3D Halftoning has been a pplied to material com¬ 
position using 3D erro r diffusion filters [Lou and Stucki 1998| 
Doubrovski et al. 2014) and 3D dispersed-dot dithering [Cho et al.| 
2003) . For color and appearance reproduction, 3D error diffusion is 
not appropriate because material assignments closer to the surface 
will have a greater infiuence on the appearance of the object than 
material assignments deeper within the object. Thus, a 3D error dif¬ 
fusion filter would have to adjust its orientation during traversal to 
account for this and maintain a consistent orientation with respect 
to the surface. An isotropic filter would produce similar artifacts 
as are observed with isotropic filters in 2D error diffusion. In con¬ 
trast, our approach of halftoning on multiple offset surfaces within 
the object, in addition to the surface itself, results in a halftone that 
inherently accounts for the geometry of the surface. The relative in¬ 
fluence of voxels at different depths from the surface is calibrated 
in an offline process and built into an International Color Consor¬ 
tium (ICC) profile. Such an offline color calibration process would 
be very challenging for a 3D filter, because it would require cali¬ 
brating every possible surface orientation. 

2D Halftoning. Generations of researchers in the field of 2D 
halftoning focused on finding methodologies to optimally arrange 
printed dots for maximizing print quality (by preserving tone and 
structure and shifting quantization errors to the highest spatial fre- 
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quencies possible - see Section [3^ subject to technical limitations 
of pri nting systems (e.g. t he ability to accurately deposit isolated 
dots) |Lau and Arce 2001) . 

One category of algorithms, called point processes, allow a very 
fast computation of the halftone screens by thresholding pixels us¬ 
ing a precomputed threshold mask that is tiled over the 2D im¬ 
age. The traditional clustered-dot order dithering to create ampli¬ 
tude m odulated screens and dispersed-dot ordered dithering | |Bayer| 
1 1973) fall into this category. The latter technique was adapted to 3D 
printing |Cho et al. 2003) accounting particularly for dot placement 
limitations of binder-jetting systems. One drawback of dispersed- 
dot ordered dithering is that frequency components given by the 
screen period are visible resulting in cross-hatch pattern artifacts. 
Avoiding such artifacts in point process techniques, requires lar ge 
threshold masks - e .g. blue-noise ma sks [Mitsa and Parker 1992) or 
green-noise masks | |Lau et al. 1999) - which are heavily distorted 
if applied on surface manifolds with non-zero Gaussian curvature. 

A second category of algorithms, called neighborhood pro¬ 
cesses, considers a local pixel neighborhood for thresholding. 
Error-diffusion methods (see Appendi x) belong to this category and 
were pioneered by Floyd-Steinberg | |Floyd and Steinberg 1976) . 
Since then, the first error diffusion techniques affected by visual ar¬ 
tifacts from low-frequency structural patterns) were improved dras¬ 
tically producing visually pleasing halftones without artifacts. Mul¬ 
tiple modifications have been prop osed including edge enhance¬ 
ment |[Eschbach~'Sid~^no^T99TJ^_tone-dependent diffusion fil¬ 
ters (^tromoukhov 2001[|Li and Allebach 2004) , threshold mod- 
ulatio n [Zhou and Fang 2003) structu re preser\^ion iQiang et al.| 
|2009| or memory efficient processing [Chang and Allebach 2003) . 

A last category employs models of the human visual system 
to optimize halft one pattern iteratively |Agar and Allebach 2005 1 
|Pang et al. 2008) . Such methods can produce the highest halftone 
quality, but are in general computationally much more expensive 
than methods belonging to the other categories. Furthermore, these 
methods iterate over the whole data and require substantial modifi¬ 
cations for adapting them into a streaming architecture. Since iter¬ 
ative methods aim to preserve local structure, their quality advan¬ 
tages over other methods diminish with increasing print resolution. 

Due to the low computational effort and high quality of 2D error 
diffusion achieved with small diffusion filters, we decided to adapt 
it to 3D color printing. One prior work addresses error diffusion 
on a surface |Alexa and Kyprianidis 2015) . This approach operates 
on meshes and traverses the vertices based on the availability of 
subsequent moves or neighbors to diffuse error to. While this ap¬ 
proach could be applied to any graph structure, including voxels, 
it is not clear that it can be applied in a streaming architecture. To 
the best of our knowledge, we are the first to consider error diffu¬ 
sion halftoning in the context of both non-Euclidean domains and 
highly translucent materials. Moreover, we are the first to propose 
such a technique demonstrated to be applicable in practice to tens 
of billions of elements. 

3. BACKGROUND 

3.1 Light scattering and absorption in printing 
materials 

Given an arrangement of multiple non-fluorescent printing materi¬ 
als with similar refractive indexes within a shape S, light propaga¬ 
tion within this s hape can be described b y the steady-state radiative 
transfer equation [Chandrasekhar I960) (see supplemental material 
for details), which shows that the intensity of light traveling through 
the material is attenuated by absorption (dependent on wavelength. 



Flashlight turned off Flashlight turned on 


Fig. 2. Subsurface light scattering in photo-polymer printing materials. 

independent of direction) and is redistributed by scattering. Thus, 
a fraction of light entering the print at one location may be emitted 
from the surface at a different location due to scattering (see Fig- 
ure[^. If light travels through different materials its spectral power 
distribution is modulated by each material’s absorption coefficients 
and the path length within this material. 

This has multiple implications for arranging translucent (high 
scattering, low absorption) printing materials colored in cyan (C), 
magenta (M), yellow (Y) and white (W) for full color 3D printing. 

Highest Reflectance (white point): To maximize reflectance, not 
only the surface must be covered with white material, but also the 
space beneath. A significant fraction of light travels deep into the 
object due to the low absorption. Each non-white material with 
higher absorption placed into the light path would absorb light 
that would otherwise contribute to the albedo. There exists a min¬ 
imum thickness dy^ of a white layer to create almost maximum 
reflectance, i.e. for each d > dy,\ \\Rdyj - Rdh < ^^max, 
where Rd^ and Rd are spectral reflectances of white layers with 
thicknesses d^ and d, and Ai^max is the maximally acceptable re¬ 
flectance difference. 

Lowest Reflectance (black point): To minimize reflectance, the 
printing material arrangement must ensure that light in the whole 
visible wavelength range is maximally absorbed. This can only be 
achieved by appropriately arranging C, M, Y materials on the sur¬ 
face and beneath. The minimal layer thickness d^ of such an ar¬ 
rangement for creating almost minimal reflectance (deviating only 
by Ai^max from minimum) satisfies d^ < dyj because the color 
materials’ absorption coefficients are larger than that of white (as¬ 
suming similar scattering properties of all materials). 

Therefore, not all voxels in the shape addressable by the printer 
need to be considered. We can drastically reduce computation by 
filling all voxels within the shape by default with white (maximiz¬ 
ing reflectance of the white point) and compute the arrangement of 
remaining materials only within d^ of the surface. Although mate¬ 
rial placement with a distance d G [db,dyj] from the surface may 
impact reflectance, it affects only very bright colors which are also 
reproducible by color halftones in upper layers. While theoreti¬ 
cally one could der ive d t based on measurements of the material 
parameters. Section [6T] shows how selecting a maximum distance 
<^max < db from the surface allows an empirical trade-off between 
the number of voxels that must be addressed and the achievable 
gamut volume. 

Optical Dot Gain: Voxels occupied by a colored material and 
surrounded by white material appear larger because the spectral 
power distribution of light traveling through both materials is mod¬ 
ulated mostly by absorption of the colored material. T his phe¬ 
nome non is known in 2D printing as optical dot gain [Rogers] 
|1997[ . The more colored voxels are stacked on top of each other 
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Fig. 3. Printing with translucent materials: severe dot gain and increased 
contrast by printing more slices. Two small cyan printed squares (indicated 
by a thin dashed line), both 8 dots or ~ 340 fim wide at 600 DPI. Left: the 
square is printed 9 slices deep into surrounding white material. Right: 36 
slices deep. Also visible is the surface topography due to printing roller. 

aligned with the surface normal the more laterally scattered light is 
modulated by the colored material. Thus, the colored area on the 
surface appears not only darker, but also larger as shown in Fig¬ 
ure Algorithms that use surface halftones to assign materials to 
inner voxels must consider this phenomenon to avoid dot gain arti¬ 
facts particularly visible in bright areas by isolated disturbing large 
dots. Section [6l2] describes how independently halftoning layers al¬ 
lows to decorrelate material assignments in different layers, which 
greatly reduces dot gain artifacts in translucent materials. 

3.2 Perceptual aspects in color printing 

Each printing system is limited in reproducing reflectances due to 
the small number of available printing materials. Still, accurate 
color reproduction is possible by material arrangements creating 
physical errors for which the human visual system (HVS) is insen¬ 
sitive. Two perceptual aspects are particularly important in printing. 

Contrast perception: Contrast sensitivity functions (CSFs) de¬ 
scribe the HVS’s sensitivity to achromatic and chromatic con¬ 
trasts as a function of spatial frequency, orientation and lumi¬ 
nance |Mullen 1985||Campbell et al. 1966|[Van Nes and Boumani 
|1967) . The achromatic CSF has bandpass and chromatic CSFs have 
lowpass characteristics. Each CSF monotonically decreases for fre¬ 
quencies higher than 10 cycles/degree resulting in an apparent blur¬ 
ring of high frequency contrasts. This is used in printing by arrang¬ 
ing printing materials in high frequency patterns (halftones) to cre¬ 
ate various color shades from only a few colored printing materials. 

Color perception: Retinal processing reduces the information 
content of a spectral stimulus to on ly three values-the cone 
responses |Wyszecki and Stiles 2000] . Therefore, different re¬ 
flectance spectra viewed under the same illuminant match visu¬ 
ally if corresponding spectral stimuli yield similar cone responses. 

By agreeing on the viewing illuminant, accurate color printing is 
possible by material arrangements that reproduce cone responses 
instead of reflectances. A drawback of this metameric color re¬ 
production is a likely mismatch between prints and originals for 
other illuminants. CIE colorimetry pro vides two standa rdized color 
spaces for metameric printing (see e.g. |Fairchild 2013| ): CIEXYZ, 
which is linearly related to cone responses; and CIELAB, an oppo¬ 
nent color space derived from CIEXYZ that allows simple access 
to color attribute predictors for lightness, chroma and hue and is 
perceptually more uniform. Color control in metameric printing is 
done by separation-voXdiXmg CIEXYZ or CIELAB values to mate¬ 
rial arrangements. Since this is generally impossible for all colors 
due to gamut limitations of the printing system, a preceding gamut 
mapping method must transform all colors into the device gamut 
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aiming t o minimize per ceptual errors between original and repro¬ 
duction IMorovic 2008| . 


4. OVERVIEW 

Given the resolution of current 3D printers, and expected increases 
in both resolution and build space, we need a streaming model 
for computation, wherein only localized data and computation are 
needed. Figure shows such a streaming framework, which takes 
as input shape and color information, and produces printer com¬ 
patible data in which a single material is assigned to each voxel. 
The error diffusion and material assignment algorithms described 
in Sections|^and[^are our core contributions and take place in the 
parts colored reddish. Each aspect of our framework is described 
briefly below. 

Shape + Color: The input to our pipeline is a shape S C de- 
flned by a surface dS (e.g. tessellated mesh) with attached color in¬ 
formation f : dS ^ C (texture data or per-vertex colors), where C 
is specifled by an invertible transformation from CIEXYZ deflned 
for an illuminant (e.g. CIED50) and a set of color matching func¬ 
tions (CMFs) (e.g. CIE 1931 CMFs). Exa mples for C are standard 
RGB color spaces ISiisstrunk et al. 19^ . By default we interpret 
RGB data as sRGB. 

Voxelization: The axis-aligned bounding box of our shape B(S) 
is discretized with a flxed, regular grid of voxels B at the resolution 
given by the printer speciflcation. We slice the surface and identify 
interior, V = H 23, and surface voxels, Vo CV = dS HB, where 
a slice refers to a 2D array of voxels with a constant ^-coordinate. 
The full slice is denoted by B{s) C B for slice s with center z- 
coordinate and Va(s) = Va fl B{s) denotes the surface voxels 
in slice s. Our voxelizer generates slices in chunks of 100 slices- 

3 mm for our hardware setup-depending on the number of layers 
and layer thickness as discussed in Section until the object is 
completely sliced. Each chunk proceeds through the pipeline until 
materials have been assigned, and printer-speciflc output has been 
generated. Then the next chunk is generated. 

During the voxelization process, we assign colors to the surface 
voxels; abusing notation slightly, we redeflne f : Va ^ C as a 
function on the surface voxels. We sample the texture data using 
trilinear interpolation (mipmapping) with the level-of-detail com¬ 
puted per surface voxel. This is key to avoid aliasing-while there 
are no viewpoint-dependent effects as in rendering, the texture map 
is in general non-uniform over the surface. 

Color Management: In this stage, which includes gamut map¬ 
ping and color separation, we map color data t o printer tona l val - 
ues using ICC color management jlCC 2010| (see sectio n [7^ : 
p : C ^ T, where each element of the tonal value space T cor- 
responds to the amount of available printing materials required to 
best reproduce a desired color. By function composition we attach 
one tonal value vector to each surface voxel: g = p o f : Va ^ T. 
Layer Construction: To account for the translucency of the ma¬ 
terials, see Figures and we transfer the tonal values stored in 
the surface voxels to the interior voxels within a distance dmax of 
the surface. See Figure for an example. These voxels have the 
greatest influence on the appearance of the printed object. From 
these voxels, we extract a set of layers as isosurfaces of distance to 
the surface of the object. The number and thickness of the l ayer s 
determine dmax, and these are chosen as described in Section |0] 
Halftoning: We then treat each layer as a distinct surface with tonal 
values, and halftone them independently using the surface traversal 
algorithm described in Section and the error diffusion and mate¬ 
rial assignment algorithm described in Section 
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—> 

• push tonal values inside 
and extract layers 
(Section 6.1) 

—> 

• consistent traversal 
(Section 5) 

• error diff. + material 
assign. (Section 6.2) 


Fig. 4. Flow-chart of our processing pipeline. 


Material Arrangement: The output of the halftoning is a voxel- 
level material arrangement m \V ^ M, where the set of materials 
M is discrete and small-4 in our case. We store m per-slice, and 
convert it to a printer-specific format to send to the device. 

5. CONSISTENT SURFACE TRAVERSAL 

Classical error diffusion algorithms fo r 2D images, e.g. | |Floyd and| 
|Steinberg 1976[[Ostromoukhov 200 1| , rely on the partial ordering 
property of subsets of Euclidean domains, to define a traversal or¬ 
der, eg. raster scan or serpentine scan order, in which error is never 
pushed to pixels that have already been halftoned. 

On a surface, which has non-zero curvature and is in general not 
a topological disc, such a partial ordering property does not ex¬ 
ist without introducing a seam or a singularity, a consequence of 
the Hairy Ball Theorem. Any traversal scheme will eventually en¬ 
counter an element that has already been halftoned, requiring either 
a stop and restart, or a change in direction. Excessive stopping and 
restarting can lead to poor distribution of error, and we design a 
traversal that avoids this. 

To map anisotropic error diffusion filters on to the surface re¬ 
quires a local coordinate system of the tangent plane at a given 
point, which is fixed with the surface normal and one additional 
direction. Mapping a filter to the surface so that it is consistently 
oriented for adjacent points is equivalent to parameterizing the tan¬ 
gent space of a surface in a consistent way, or finding smooth vec¬ 
tor fields on a surface, which is in general an ill-posed problem and 
finding an optimal solution requires considering the entire surface. 

Instead, we propose a voxel-level traversal algorithm that main¬ 
tains a consistent orientation of the filter, traverses the surface vox¬ 
els in long runs allowing error to accumulate and diffuse evenly, 
and does so using only local information. Our traversal avoids situ¬ 
ations where voxels only distribute or obtain error during diffusion. 

We maintain a consistent orientation of the filter by traversing the 
voxels in a consistent direction, and using the traversal direction to 
fix a coordinate system in the tangent plane. While our traversal 
provides no theoretical guarantees, it performs well in practice. 


5.1 Geometry of Voxel Slices 

As a consistent direction, we always travel clockwise or counter¬ 
clockwise about a line parallel with the positive z-axis. For brevity, 
we describe only the counterclockwise case. Note that we want to 
traverse counterclockwise, not about a global axis, but rather about 
a line through the center of the local connected component of vox¬ 
els. Depending on the geometry of S, when it intersects a slice 
plane ^ and is discretized into a voxel representation, the re¬ 
sulting voxels in slice s may belong to multiple connected compo¬ 
nents. An example of a slice with multiple connected components 
can be seen in Figure 

If the surface normal is close to vertical (almost parallel with 
the positive or negative z-axis), then Va(s) will be multiple voxels 
wide. These surface voxels Va(s) may enclose interior voxels in 
slice s, and will be themselves enclosed by exterior or empty voxels 
in slice s. Figure|^shows a 2D cross-section of this. 

If we consider a part of the shape, which forms a connected com¬ 
ponent in consecutive slices, and where the surface faces almost in 
the negative z direction, but slopes or curves slightly upward, we 
have the following situation, which is shown in 2D in Figure 
When projected into the xy-plane, the surface voxels of Va(s) will 
be located inside the surface voxels of Va (s+1). If we consider two 
surface voxels within the same connected component in slice s, the 
innermost should be traversed first, since they are next to Va(s — 1) 
and have already received error diffused from slice s — 1. The out¬ 
ermost should be traversed last, which are next to Va(s + 1). The 
situation is reversed for upward facing surfaces. 

5.2 Traversal Algorithm 

We construct an undirected graph for traversal as follows: for each 
surface voxel v G Va, we connect it with an edge to each other sur¬ 
face voxel within the 3x3x3 voxel window W(v) centered on v. 
We denote this set of one-ring neighbors M{\) C Vo- We traverse 
the voxel representation of the surface slice by slice, analogous to 
halftoning an image row by row, by selecting the next voxel from 
the set of neighbors that are within the same slice, A4(v) C A/’(v). 


ACM Transactions on Graphics, Vol. 35, No. 1, Article 4, Publication date: December 2015. 






























6 


A. Brunton et al. 




Fig. 5. Geometry of surface voxels in consecutive slices for a down-facing 
surface. 


Error is always diffused upwards to the next slice, and to other sur¬ 
face voxels within the slice, which have not yet been assigned ma¬ 
terials, as shown for a very small example (1 cm) in Figurej^ 

One might also consider standard graph traversals, such as 
breadth-first or depth-first, also constrained to neighbors within the 
same slice, but these do not maintain a consistent direction of travel 
and therefore only permit symmetric error diffusion filters (e.g. uni¬ 
form, Gaussian, etc.). _ 

As described in Section the traversal winds up the surface 
on a slice-by-slice basis counterclockwise about a line through the 
center of each connected component. However, finding the cen¬ 
ter of the current connected component would require a search 
over all voxels in the slice to determine the connected compo¬ 
nents, find their centers, and to which connected component each 
surface voxel belongs. Such a search would prove too expensive- 
essentially involving a pre-traversal. 

We can determine whether a potential next voxel in proper traver¬ 
sal order winds about the connected component in the correct di¬ 
rection using only two pieces of local information. The first is the 
surface normal. The second is the in-slice distance to the exterior 
of S, or distance-to-empty, denoted by 0s (v) for voxel v G Va(s) 


0s (v) 


min llv — u||i , 

ueB{s)\v{s) ^ 


( 1 ) 


which can be pre-computed for each slice-see Figure]^ for an ex¬ 
ample. We use the Li-norm for simplicity and efficiency, but an¬ 
other norm would also work. The gradient of 0s (v) always points 
to the interior of the connected component to which v belongs. This 
is important for situations where the current slice may have multi¬ 
ple disconnected components, as in Figure 

This allows us to efficiently test if a potential next voxel u G 
A/’s (v) continues in the same direction around the connected com¬ 
ponent by computing (u — v) x V0s(v). For example, if we are 
traversing counterclockwise, we can rule out all neighbors of v for 
which this cross-product has a negative ^-component. We denote 
the subset of A/’s (v) satisfying this direction criterion A/^ (v). 

Referring to Figure and recalling the observations of Section 
113 we get the following additional traversal criterion. If the sur¬ 
face faces downward, we wish to stay as inside-SiWSLy from exterior 
voxels-as possible. If the surface faces upward, we wish to stay as 
outside as possible. Hence we choose the neighbor 


w = arg max 0^ (u) (2) 

ueA/'x (v) 

for down-facing surfaces, and 

w = arg min 0^ (u) (3) 

ueA/'x (v) 


Otherwise. 


Fig. 7. A visualization of traversal order for a slice containing both upward 
and downward facing parts: (a) slice plane intersecting the object; (b) tonal 
values at the surface in that slice; (c) color-coded traversal order with color¬ 
coding legend (bottom=first, top=last). 


Because the surface orientation is determined locally using the 
surface normal, the traversal adjusts as it traverses a slice with both 
up- and down-facing segments. Figureshows the traversal order 
for such a slice. 

5.3 Boundary Cases 

Start point selection. We also use the distance-to-empty, in combi¬ 
nation with one other local quantity, to select the start point in each 
slice. First, potential starting voxels (voxels in the slice, which have 
not yet been traversed) are filtered based on their error diffusion 
count, or number of times they have had error diffused to them from 
previous slices. We select the subset of voxels with the largest er¬ 
ror diffusion count, and subsequently select from this subset using 
the distance-to-empty. If the ^-component of the surface normal is 
negative, the start point is the point, which maximizes the distance- 
to-empty. If the ^-component of the surface normal is positive, the 
start point is the point, which minimizes the distance-to-empty (i.e. 
one of outer-most points). 

Once the start point and traversal direction are selected, the direc¬ 
tion is tested to ensure it is not a dead end. If no voxel is available in 
that direction, the direction is reversed. This reduces the occurrence 
of single voxels with few neighbors to push error. 

Birth and death of components. During the slicing process, new 
parts of the shape are encountered, creating new connected compo¬ 
nents of voxels in a given slice. When a new component is encoun¬ 
tered in a given slice, the surface voxels will form a disc in that 
slice, meaning it can be traversed like an image. We handle these 
cases specifically as follows. Rather than spiral out, as with other 
downward facing surfaces, we select a global axis, x or y, and per¬ 
form a 2D serpentine scan in this direction. This avoids situations 
where a long strip of voxels mostly has error pushed away from it, 
creating a start-up artifact noticeable for low tonal values. We also 
use a 2D serpentine traversal for components, which are ending in 
the current slice, rather than spiraling in. 

We detect the birth of a component when all possible starting 
points have an error diffusion count of 0, and the death of a com¬ 
ponent when all possible starting points have no neighbors in slice 
s + 1. 

Serpentine Surface Scan. When we can no longer traverse in the 
same direction locally (all u G A/’x (v) have been traversed), we re¬ 
verse direction, as in the serpentine traversal scheme from 2D error 
diffusion, which is shown to provide significant improvement over 
raster-scan traversal in images |Ulichney I987[ [Chang and Alle-| 
[bach 200^ . Further, when selecting a new start point, we set the 
traversal direction opposite the traversal direction, in which error 
was last diffused to the new starting voxel. 
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Fig. 6. Traversal of a surface maintaining consistent orientation of an error diffusion filter. Left: a rendering of a 1 cm voxelization of a head model before 
and after halftoning. Middle: a close-up of the halftoning process up to a given slice. Right: a cut-out from the middle of how a 3-component filter moves as 
each voxel is quantized. Note: the halftone is shown as opaque and diffuse for illustration purposes, and does not reflect the printed appearance. 


5.4 Mapping the filter 

Once we have chosen the next voxel to traverse, we can establish a 
local coordinate system at the current surface voxel, with which we 
can align the error diffusion filter and the neighbors of the current 
voxel. As in image error diffusion, one non-zero filter element is 
aligned with the next voxel to be traversed. This gives us a direc¬ 
tion, which together with the surface normal, gives us the necessary 
coordinate system for the tangent plane at each surface voxel. De¬ 
pending on whether the direction is clockwise or counterclockwise, 
we use either a left-handed or right-handed coordinate system. Fig¬ 
ure!^ shows how an error diffusion filter with 3 non-zero elements 
is mapped to the neighbors in the tangent space of a voxel as it 
traverses the surface voxels. 

This is closely related to decal mapping methods such as discrete 
exponential maps | |Schmidt et al. 2006) , and to parallel transport on 
manifolds, and for large decals exponential maps provide a more 
general solution. However, the filters used for error diffusion are 
typically small enough that the neighbors can simply be projected 
into the tangent plane while preserving distances. 

We let the error diffusion filter live in the local tangent coordi¬ 
nate system, orthogonally project the neighbors onto the tangent 
space, and find the best alignment of neighbors and non-zero filter 
elements. We align neighbors with filter elements using symmetric 
closest points to avoid a single neighbor being assigned to multi¬ 
ple filter elements. That is, both the filter element and the neighbor 
must be closest to each other. 

5.5 Analysis 

The distance to empty can be computed in time linea r in the num- 
ber of voxels in each slice using a distance transform | |Felzenzwalb| 
land Huttenlocher 2004) . While this means over the whole print job, 
every voxel in B must be visited, the constant hidden in the 0( ) 
notation is very small and the distance transform is highly efficient. 
Selection of the starting point requires a search over all unprocessed 
surface voxels, but for non-pathological cases this will only be nec¬ 
essary a small number of times per slice. Each traversal operation 
requires a constant number of operations, as the number of neigh¬ 
bors is constant. Normals can be computed in 0(| Va(s) |) per slice 
using a signed dis tance field as an implicit surface definition, as de¬ 
scribed in Section [O] Thus, the total time complexity to fully pro¬ 
cess S is linear in the total number of voxels 0{\B\) with a small 
constant, which is very well suited to a streaming architecture. 


6. LAYERED HALFTONING 

Applying a halftone only to the surface voxels will result in a very 
low contrast print, due to the translucency of the materials. We can 
see this from Figure which shows how increasing the depth of 
color material assignment increases the contrast. In Figure we 
see how increasing the depth of color assignment increases the 
color gamut. For this reason, we introduce layered halftoning, in 
which color data is transferred from the surface to additional iso¬ 
surfaces, or layers, within the object, which are then independently 
halftoned usin g the traversal scheme presented in Sectionand an 
adaptive fi lter ||Ostromoukhov 2001) and threshold error diffusion 
algorithm |Zhou and Fang 2003f 

6.1 Layer Construction 

Layered halftoning begins by transferring the surface tonal values g 
from Vd to a subset of V within a distance timax of the surface. All 
interior voxels more distant from the surface are assigned a white 
material to maximize reflectance as explained in Section This 
requires both many distance computations and a transfer operator 
T : Va X T ^ V x T, which maps functions on Va to functions 
on V. 

To compute the distance field d : B ^ of voxels to the 
surface, we leverage the fact that we are only interested in voxels 
within a given distance c/max- This allows us to construct a 3D mask 
containing these distances at the printer resolution. The distance 
field d{\) is initialized with a large value > o^max- Distance 
computation then proceeds by moving this mask over the surface 
voxels Va, and writing the mask value to the voxel at the corre¬ 
sponding offset, if that value is less than the value already there. 
While this approach is efficient only for relatively small distances, 
it has the advantage of being embarrassingly parallel and requiring 
only a less-than operation at each voxel. 

As a transfer operator, we simply take the tonal value g(u) of 
the nearest surface voxel u, allowing the transfer operator to piggy¬ 
back on the distance computation at minimal extra cost, i.e. interior 
tonal values are assigned as 

g(v) = g(u) such that u = argmin ||v — sjL (4) 

seVa 

for V G V. Aside from being easy to compute, this has the benefit 
of preserving high-frequencies in g. Figure shows the process of 
assigning sRGB color values to surface voxels, converting them 
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Fig. 8. Layered halftoning process: (a) slice plane 2 ; = 2 :^ intersecting the object; (b) the sRGB values in Va(s) converted to tonal values g via the map p; 
(c) transferred tonal values g inside the object; and (d) material assignments m (black indicates no material). 


to tonal values, transferring the tonal values to interior voxels and 
applying error diffusion to get the final material assignment. Note 
that because tonal values are pushed to the interior in 3D, interior 
tonal values may include surface tonal values from different slices. 

Although there are inevitably discontinuities in g at the cut loci 
of the distance field d (voxels v where multiple surface voxels are 
equidistant), we found this not to be a problem as the tonal values 
mostly vary smoothly without significant correlation between high 
frequencies in g and high curvature of the surface. In this case we 
simply assign tonals in a first-come fist-served manner. Similarly, 
although areas with varying curvature result in slight tone shifts 
from g to g, in particular for areas of high positive mean curvature, 
we did not find this to be a problem, even for small prints. 

We extract layers of voxels with tonal values from V by defin¬ 
ing isosurfaces of d with which to segment V. As isosurfaces we 
choose dj^ — ir, where r is the layer thickness, which we choose to 
be the voxel size along the dimension of minimum resolution. We 
use ^ = — 1, where L is an integer co nstan t, which de¬ 

fines (imax = Lr. Ideally timax = (see SectionjTTJ to minimize 
reflectance for maximizing contrast. However, we make a tradeoff 
between computational effort and gamut volume and set L = 12 
for all prints shown in this paper. Figurej^shows color gamuts com¬ 
puted using L = 3,6,12,18 with the gamut volume as a fraction 
of the volume of sRGB. We can see how the volume increases with 
the number of layers. However, we found that the black point of the 
18-layer profile, L* = 32.14, is close to the minimum black point 
achievable with current CMYW materials-L* = 31.55 for a 3 cm 
cube of full CMY mixture. Thus, little additional gamut is likely to 
be gained by more layers. 

Voxel layer i is then defined as the set of voxels between isosur¬ 
faces i and ^ + 1, 


= {v G V : di < d{\) and 3 u G W(v) : d{u) < di} (5) 


where the second condition is to ensure that the layers form thin 
voxel approximations of isosurfaces as muc h as possible, and W 
is the voxel window defined in Section 15.21 Due to different res¬ 
olutions along each axis, — di may be multiple voxels thick 
depending on the orientation. In the case of Vo, this condition is 
replaced by 3 u G yV(v) such that u ^ V. 

Note that we can very efficiently compute normals for all isosur- 
face s, required for traversal of the error diffusion filter (see Section 
5.2 1 : n(v) = V(i(v), v G V^, where 


d(u) = 


—(i(u) if u G V 
d{u) otherwise 


( 6 ) 


is a signed distance field, for u e B. We use finite differences to 
compute Vd. 


6.2 Error Diffusion and Material Assignment 

We treat each voxel layer as a set of surface voxels, traverse them 
as described in Section[^ and halftone them as follows. Each tonal 
channel is halftoned independently using a stan dard error diffusion 
filter, mapped onto V^ as described in Section |5.4| and in the Ap¬ 
pendix. At a given voxel v G V^, this results in a halftone vector 
h(v) G {0,1}^ where T = dim(T). Thus, multiple materials may 
have a tonal value of 1 indicating they should be assigned to v. 
Since we can only assign a single material m{\) to each voxel, we 
apply the following tie-breaking scheme. 

For each tonal channel t G {1,..., T}, we maintain a tie-breaker 
error, c(t), which is initialized to c(t) = 0 for all t at the start 
of each slice. If at a given voxel, one or more tonal values are 
halftoned to 1, the tonal channel t* = argmaxtc(t) with the 
largest tie-breaker error is declared the “winner” and we assign the 
corresponding material to the voxel. We then reset c{t*) = 0 and 
increment c{t) for all other t. If there is no t such that h(v) [t] = 1, 
we assign white material to the voxel. 

Note that this tie-breaking scheme is independent of the per-tonal 
error diffusion, and the error diffusion does not know when a tonal 
is quantized to 1, but the material is not assigned. This in general 
leads to a slight loss of tone, but is mitigated by the large number of 
voxels, and it is any way a ccounted for by the color management, as 
described in Sectionbecause we characterize the printer using 
targets printed with the same algorithm. 

Note that errors are only diffused within the same layer, and each 
layer is halftoned independently. This has the following benefits. 
First, it allows different layers to be halftoned in parallel. Second, 
it further avoids the issue of resampling a halftoned signal, which 
would be necessary if the halftone was performed on the surface 
and then transfered to the interior. This is problematic because the 
halftoned signal is by design high-frequent and resampling it in¬ 
troduces artifacts. Third, we can use the thres hold modulation of 
the halftoning algorithm |Zhou and Fang 200^ to avoid inter- laye r 
correlations, which can create dot gain artifacts (see Section [3T] ). 
We do this by seeding the pseudo-random number generators used 
for the threshold modulatio n diff erently for each slice and layer. 

As discussed in Section |6.1| due to differences in resolution 
along the different axes, some interior voxels will fall in between 
layers, and not be directly halftoned. For these voxels, we simply 
assign the material of the nearest voxel that does belong to a layer. 


7. EXPERIMENTS 
7.1 Printing Setup 

We use an Objet500 ConnexS from Stratasys |Stratasys 2014) , for 
which the vendor has provided an interface allowing us to specify 
material assignments at the printer resolution. As the ConnexS al¬ 
lows only 3 model materials and 1 support material, we expanded 
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b* a* b* a* b* a* b* a* 


0.269 0.322 0.372 0.409 

Fig. 9. Increasing number of layers results in increasing gamut volume. From left to right: 3, 6, 12 and 18 layers. The gamut volume is indicated below each 
image as a fraction of the volume of sRGB gamut shown transparent. 


the possible gamut by coloring the support material with a yellow 
acid dye. Thus, we print with white, cyan, magenta and yellow ma¬ 
terials. Care must be taken that there are no large agglomerations 
of support within the object, which would weaken the structural 
integrity of the print. Fortunately, as our halftone is by design high- 
frequency, this is avoided by scaling the yellow tonal value by a 
constant factor between 0 and 1. We found 0.3 to provide both 
a chromatic yellow and structurally stable prints. Our layer-based 
approach also allows us to prevent support from being present in 
the outermost layer (the surface), by applying a layer and tonal de¬ 
pendent soft-thresholding to the tonal values. For error diffusion 
we use a tone-adaptiv e filter ||Ostromoukhov 20'0T| and threshold 
modulation technique [Zhou and Fang 200^ 

7.2 Color Management 

Specifying the color-to-tonal transformation p requires at least two 
ICC profiles: An input profile that specifies texture-data colori- 
metrically and a printer profile that includes gamut mapping and 
color separation. Generating a printer profile involves a colori¬ 
metric characterization of the printer, i.e. we need to predict the 
printout’s CIEXYZ color for each tonal value vector in 'T . For 
this, the tonal value space is sampled and the color of correspond¬ 
ing printouts is measured and used to fit a colorimetric printer 
model. We used a fully empirical model, which interpolates be¬ 
tween color-measured printouts of densely sampled tonal values- 
for our CMYW printer we use a uniform 8x8x8 sampling of 
the channel-wise linearized CMY tonal value cube. Linearization 
is perform ed employing the broad band Murray-Davies model ac¬ 
cording to fWyble and Berns 2000| . This target is printed using the 
algorithm described in Sections[^nd[^exactly the same algorithm 
as used to generate the results shown in Figures [TT]and[^ without 
color management as the input are direct tonal values. A picture of 
the resulting printout is shown in Figure(left). 

Measuring colors of such highly translucent printing materials is 
a challenging problem. It can be shown that measurements made 
by spectrophotometers used in the graphic arts community are sys¬ 
tematically biased towards lower reflectance due to subsurface light 
transport away from the detection area. For our color measure¬ 
ments, we used a bidirectional 0°/45° measurement geomet ry, an 
almost colorimetric DSLR camera (Vora value of 0.9455 |Vora| 
land Trussell 1993] ), diffuse broadband illumination simulating 
CIED50, flat fielding (to account for optical path variations be¬ 
tween light source, printed color patches and camera), and a poly¬ 
nomial approach to map camera RGB values to CIEXYZ. Color er¬ 
rors determined by spectroradiometric measurements are within the 
interinstrument-variability of spectrophotometers used in graphic 
art s for opaque mater ials. Eor more details, the reader is referred 
to [Arikan et al. 2015| . 

We evaluated the effectiveness of color management as follows. 
We took a standard color checker and measured each of the 24 




Fig. 10. Left: color characterization target with 12 layers, which was used 
to compute the ICC profile use to print the results in this paper. Right: pre¬ 
dicted and measured color patches for a color checker; the left of each patch 
is the predicted color (gamut mapped) and the right is the measured color 
of the printed patch. 

patches with a spectrophotometer. We then mapped each of these 
colors into our printer gamut (12-layers, Figure second from 
right), and printed planar patches. We then measured the printed 
patches with the almost colorimetric DSLR cam era, and computed 
the CIEDE2000 [CIE Publication No. 142 20011 differences to the 
gamut mapped colors-i.e. the colors predicted by our empirical 
printer model. We observed a median error of 2.2, a mean error 
of 2.3, and minimum and maximum errors of 0.5 and 6.1. The pre¬ 
dicted and printed colors are shown in Figure(right). Note that, 
even though the proposed half-toning algorithm was used to print 
the planar patches, this evaluates only the empirica l pri nte r mo del 
and not the halftone, which is evaluated in Sections [73] and lTT] 

7.3 Evaluation 

Performance. Our non-optimized implementation supports multi¬ 
threading for some computations (for example the transfer operator 
and halftoning each layer in parallel), but does not exploit any GPU 
capabilities, although several aspects would lend themselves well to 
GPU implementation. Table [T| shows performance data for several 
prints: the number of voxels, both those assigned a material (|V|) 
and in the build volume (|i3|), computation times include the time 
to deliver the first chunk of slices to the printer and the total, and 
the print times. The results were collected on a standard desktop 
PC with an Intel Core 17-4770 processor and 32 GB of memory. 
Tone Preservation. It is in general difficult to evaluate how well a 
3D print reproduces the color of the input texture, due to the large 
change in appearance resulting from different illuminants and illu¬ 
mination distribution in co mbin ation with the curved surface. As 
discussed further in Section [T6| the printing materials we use are 
highly color inconstant, meaning the perceived color changes un¬ 
der different illuminants. Further, typically the textured model is 
not rendered using the reflectance properties of the printing mate¬ 
rials. In the case of scanned objects, evaluating with respect to the 
original object is meaningless, since the capture systems are typi¬ 
cally not color calibrated. 

Therefore, we evaluate the quality of tone preservation, rather 
than the similarity of perceived colors, as follows. We compute the 
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Model 


Bald head (20 cm) 

Figurefl^ 

Ruslan (lo cm; 

Figurefl^ 

Kuslan (b cm; 


FigurefT^ 


FigureJT^ 

Apple cm; 

Figurefl^ 


Table L Performance details of our method. 


Voxel Count 
|V| \B\ 


Compute Time 
first chunk total 


Print Time 


Material-Tonal RMSE 
M Y 


8.6 X 109 2.49 X IQi 


1 min 


6 h 10 min 


21 h 


0.0038 0.0074 0.0074 


5.7 X 109 1.36 X IQio 


1 min 


4 h 18 min 


11 h* 


0.0097 0.0038 0.0107 


2.12 X 10® 5.07 X 10® 


16 min 


2h 


0.0094 0.0041 0.0120 


5.0 X 109 1.67 X 10i° 


1 min 


4 h 55 min 


20 h 


0.0056 0.0041 0.0080 


1.8 X 10^ 3.7 X 10^ 


30 s 


1 h 


10 h 


0.0056 0.0030 0.0308 


See text for a description of how we compute the tonal errors. *Print job terminated after 69% of slices. 


W 

0.0184 


0.0236 


0.0243 


0.0173 


0.0247 


root mean squared error (RMSE) between material assignments 
m and tonal values g, which are given for several models in the 
last column of Table |I| First, we compute material fractions for 
each slice-the number of voxels assigned a given material in the 
slice divided by the total number of voxels v G V(s) such that 
d{\) < drnax- Thcu wc determine expected material frac tions from 
the tonal values g(V(s)) using the Demichel equations |Demichel 
[T^ . For example, we compute the expected material fractions 
for cyan and white from tonal values as 

Ec = C(1 - M)(l - y) + |CM(1 - Y) 


+ \C{1- M)Y +\CMY 

(7) 

II 

1 

o 

1 

S 

77 

1 

(8) 


respectively, where C, M and Y are the average tonal values over 
the same set of voxels. The formulas are similar for magenta and 
yellow. Finally, we compute the errors as the difference between 
the expected material fraction and the actual material fraction for 
each tonal channel. We see that our algorithm preserves tone very 
well, with most errors being less than 1% of the tonal value range. 

7.4 Visual Experiment 

Visual experiments were conducted to assess the structural quality 
of the proposed halftone method, in particular the level of graini¬ 
ness and the visibility of structural artifacts. For this, we generated 
a test surface with the following function 

z{x, y) — (0.5 cos(3r) + 0.5) (9) 

where r — a — 2i7/3, = 5 cm with both x and 

y ranging from —5 cm to 5 cm. We then applied two textures to 
this surface, as shown in Figure pT] one with smooth color gradi¬ 
ents and one with high-frequency patches; the colors were selected 
to uniformly cover the a*b*-plane of lightness L*=70 and gamut- 
mapped using the absolute-colorimetric intent of our ICC profile. 

The printed surfaces were placed on a colorimetrically charac¬ 
terized display next to their rendered versions (gap of approx. 1.25 
cm - the left/right position of rendering and print was randomly 
chosen for each subject). We used an Eizo ColorEdge cgSOlw 30- 
inch display with a resolution of 2560 x 1600 pixels. The printed 
surfaces were 10 cm x 10 cm, which covers 400 x 400 pixels 
on the screen. The screen was placed horizontally and illuminated 
from directly above using a JUST NORMLICHT LED Color View¬ 
ing Light M viewing booth at a distance of 1 m from the screen 
to provide a diffuse illumination. The illuminant used was LED- 
simulated CIED50. A chin-rest was placed so to obtain 60 pixels- 
per-degree (ppd) at the screen center (85 cm distance) and a 45° 


viewing angle to the screen. The renderings were done to approxi¬ 
mately match this viewing condition and the luminance level of the 
print. Figure [TT] shows pictures taken with an almost colorimetric 
DSLR camera of what the subjects viewed. A detailed image and 
diagram of the setup can be found in the supplemental material. 

30 naive subjects with normal visual acuity - tested with the 
Snellen test - participated in the experiments (10 female and 20 
male; average age of 30 years). For each of the two textures, a total 
of 10 distorted textures were created by adding zero-mean Gaus¬ 
sian pixel-noise with standard deviations of 1,..., 10 CIEDE2000 
units, simulating different levels of graininess. In a first experiment, 
subjects were asked to select the rendering with the distorted tex¬ 
ture matching the graininess of the print. In a second experiment, 
they had to judge the level of structural artifacts in the print - show¬ 
ing the noise-free rendering as reference - on a quality scale of 0-5: 
0 (not visible), 1 (visible but not disturbing), 2 (slightly disturbing), 
3 (disturbing), 4 (very disturbing), 5 (extremely disturbing). 

Results of the experiments are shown in Table For both 
textures, the average perceived graininess level corresponds to a 
noise standard deviation smaller or equal to 2 CIEDE2000 units. 
To put this value into perspective, we computed the pixel-wise 
CIEDE2000 differences between the noise-free and distort ed tex- 
ture (noise std. of 2 CIEDE2000 units) using S-CIELA B [Zhang 
and Wandell 1996| with a visual resolution of 60 ppd (Johnson 
and Fairchild 2003| . The 99th percentile of CIEDE2000 errors is 



Fig. 11. Pictures of the test surfaces taken from the chin-rest position in 
our experiment; printed (left) and rendered (right). Moire on the renderings 
is caused by the acquisition process. 
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Table IL Results of the visual Experiments. 
Experiment I (graininess level [noise std. in CIEDE2000]) 



mean 

std 

min 

max 

Smooth texture 

0.9667 

0.9994 

0 

4 

Patch texture 

2.0000 

1.2594 

0 

4 

Experiment II (structural artifacts [quality scale 0-5]) 


mean 

std 

min 

max 

Smooth texture 

1.3000 

0.7497 

0 

3 

Patch texture 

1.3667 

0.8899 

0 

4 


1.644 (smooth texture) and 1.694 (patch texture). Since this is only 
slightly above the just noticeable distance, we can conclude that 
the perceived graininess of the halftone is very low. The second 
experiment revealed that also structural artifacts do not adversely 
affect the perceived quality of the halftone: for both textures, 
average quality scores are smaller than 1.4, i.e. artifacts are visible 
but not even slightly disturbing. Note that there are some drying 
related artifacts, which are likely considered in the quality scores 
but independent of the halftone. 


7.5 Qualitative Results 

Figure shows some qualitative results demonstrating the detail 
and realism produced by our software using our hard ware and ma¬ 
terial setup. Part (a ) shows (left-to-right) a head scan | |ScanLab and| 
[TurboSquid 2013| with original texture rendered with diffuse shad¬ 
ing, the same model with the texture gamut-mapped using our pro¬ 
file, and the resulting 3D print (15 cm tall). Note how well both 
smooth tones and high-frequency details are reproduced. Part (b) 
shows the same model printed 5 cm tall to show the effects of scale 
on the prints-note how similar the features of the two prints are 
despite the di fference in scale. Part (c) shows a realistic-looking 
printed apple [TurboSquid 2010| . Parts (d) and (e) demonstrate how 
well our halftoning algorithms reproduce high-frequency details in 
the textures. Prints (c)-(e) were generated with an error in the ICC 
profile, so the colors are off in abso lute terms when compared to the 
textures. In part (d) | |Ten24 2013| , the back of the ear is rendered 
next to the printed ear, which is about 3.5 cm high. The images of 
the prints in parts (a) and (e) were taken under simulated CIED50, 
to which we color-characterized the camera. We transformed the 
raw R GB images to CIEXYZ values as described in [Arikan et al.| 
|2015[ Sec. 3.2], and transformed these values to sRGB for display. 
All other photographs of prints were taken using automatic camera 
settings. 

7.6 Limitations 

Our setup is limited primarily from the hardware side. With four 
materials we are forced to print without a black. This means 
we have difficulty reproducing dark colors and a lack of color 
constancy-a change of illuminant alters the perceived color. We 
measured the color inconstancy index (CII) between CIED50 and 
CIEFll of cyan, magenta, yellow and w hite, to be 4.65, 3.24, 5.28 
and 0 .59, respectively, using CAT02 [CIE Publication No. 159| 
|2004) for chromatic a daptation and CIEDE2000 for color differ¬ 
ence [Fairchild 201 3| . We also have no control over the translu- 
cency, as the materials have very similar scattering and transmis¬ 
sion properties. We can only use their translucency to create real¬ 
istic looking objects, not control it to a desired level. The translu¬ 
cency of the materials also results in blurring of some details; an 
interesting approach has bee n proposed to addres s this in a way 
similar to unsharp masking [Cignoni et al. 2008) . The introduc¬ 


tion of printers with more materials, and new, more opaque mate¬ 
rials would address these limitations, and open the possibility of 
combining color with desired scattering and reflective properties. 
A compelling open question is whether such complex appearance 
properties can be achieved using error diffusion-style techniques. 
A further limitation is due to the colored support material. Dur¬ 
ing the uv-curing process, support material adjacent to model ma¬ 
terial mixes and binds with the model material. This results in a 
yellow tinge on some sloping and vertical surfaces. This could be 
overcome without additional materials by a multi-pass printing, in 
which the support material is cured separately from the other ma¬ 
terials. 

Algorithmically, our approach has the same limitations as 2D er¬ 
ror diffusion methods: for extremely low tonal values, we encounter 
mild start-up artifacts-materials are not placed entirely uniformly. 
We have found this only in extreme cases. 

The curvature of the surface affects the transfer of tonal values to 
the interior layers when the mean absolute curvature H approaches 
l/d^, where is the depth of highest reflectance as discussed 
in Section equivalently as the radius of curvature approaches 
d^. High negative mean curvature (convex) regions may darken 
and high-frequency texture features that coincide with high posi¬ 
tive mean curvature (concave) regions will “bloom”. We have not 
found this to be a major problem, although such effects show up 
strongest in small prints, e.g. the nose of the small Ruslan print in 
Figure [^(b). We plan to address curvature-dependent effects in 
future work, e.g. by adjusting the tonal values according to curva¬ 
ture, which could be estimated using the signed distance field d in 
([§. The surfaces used in the visual experiment also demonstrate the 
extent to which curvature influences the perceived color. These sur¬ 
faces, shown in Figure[^ contain both concave and convex regions 
with high curvature, and relatively flat regions on the periphery. 

8. CONCLUSIONS 

In this paper, we have presented algorithms for precise and effi¬ 
cient control of material placement in multi-jet 3D printers for the 
purposes of halftoning-a critical ingredient in accurate color repro¬ 
duction. We have presented a layered halftone, which accounts for 
the high-translucency of currently available color materials, and a 
traversal algorithm for voxel representations of surfaces, which al¬ 
lows us to map 2D error diffusion filters onto a surface in a consis¬ 
tently oriented way. This allows us to leverage decades of knowl¬ 
edge from 2D error diffusion research. Our algorithms fit within an 
efficient streaming architecture, preserve tone and do not produce 
artifacts. We have further shown how our algorithms seamlessly al¬ 
low the printer’s support material to be colored, thereby expanding 
the printer’s color gamut while preserving structural stability. The 
3D color prints generated with our setup exhibit a high level of de¬ 
tail and realism. 

Due to the ability to arbitrarily configure materials with different 
optical properties through an object, multi-jet printing has tremen¬ 
dous potential for graphical 3D printing in terms of reproducing 
complex appearance properties. By providing the first full color ca¬ 
pabilities for these printers, we have provided a critical building 
block for graphical 3D printing. 

The introduction of more opaque color materials will allow our 
algorithms to print with larger color gamuts using fewer layers, re¬ 
sulting in a performance boost. Better exploiting parallelism in our 
computations is a key direction for future work as printer resolu¬ 
tions and build volumes increase. As the number of printing ma¬ 
terials increase, it will be important to consider printer models in 
characterizing the printers and computing profiles. 
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Fig. 12. Various 3D color prints generated using our software. Note that color deviations between rendering and print may have various reasons such as 
gamut limitations, goniochromatism and color inconstancy of printing materials, light field and illuminant differences between real scene and rendering, etc. 
The comparison aims to show how well texture details are preserved by our halftoning approach and that no artifacts are introduced. See the text for discussion 
of the individual parts. 
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Appendix 

Our halftoning algorithms we propose are based on error diffusion, 
in which after a pixel or voxel is quantized, the resulting error is 
distributed to nearby pixels or voxels using an error diffusion filter. 
Thus, if a pixel with a value of 0.5 is quantized to 0, the error 
is distributed to its neighbors, making it more likely that they are 
quantized to 1. Formally, we are given an input tonal signal g and an 
error signal a, which is initialized a(v) = 0 V v. The error diffused 
signal is defined as 


g(v) = g(v)+a(v), (10) 


h(v) 


1 if g(v) > t 
0 otherwise 


( 11 ) 


the resulting error is then diffused to the neighbors 

a(u) =a(u) -«;v.u(h(v) -g(v)) (12) 

where u G ^(v) is a neighbor of v, is an element of the error 
diffusion filter diffusing error from v to u. 
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9. SUPPLEMENTAL MATERIAL 
9.1 Radiative transfer 

Given an arrangement of multiple non-fluorescent printing materi¬ 
als with similar refractive indexes within a shape S, light propaga¬ 
tion within this shape can be described by the steady-state radiative 
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transfer equation |Chandrasekhar I960) 

^^^^^ = -axix)h{x,Q) (13) 

- s>,{x) (^hix, n) - £ hix, n')p^in', n)dn') 

where Ix{x,Q) is the radiant intensity at location x e S prop¬ 
agating in direction Q for wavelength A G [380, 730] nm (vis¬ 
ible wavelength range), ax{x) > 0 is the spectral absorption 
coefficient, sa(x) > 0 is the spectral scattering coefficient, 
is the sphere and px{^',^) is the scattering function satisfying 
f^)dQ'dQ = 1. Solving the radiative transfer equa¬ 
tion with an appropriate boundary condition at the surface gives us 
the radiation that is diffusely emitted from the shape. Adding the 
radiation directly refiected from the surface due to Fresnel reflec¬ 
tion, gives us the total radiation emitted from the object’s surface. 

Eq. {13} shows that the intensity of radiant energy traveling 
through the material is attenuated by absorption (note that ax{x) 
depends on the wavelength but is independent of the traveling di¬ 
rection) and is redistributed by scattering (scaled by sa(x)). Thus, 
a fraction of light entering the print at one location may be emitted 
from the surface at a different location due to scattering. If light 
travels through different materials its spectral power distribution is 
modulated by each material’s absorption coefficients and the path 
length within this material. 

9.2 Gamut Volume Visualizations 

Figure[^ shows the gamut volumes from the paper projected onto 
the a*b* plane to provide another perspective on how the gamut 
volume changes as layers are added. 

9.3 Visual Experiment Setup 

Figure shows the physical setup we used to conduct the visual 
experiment. Subject placed their chin on the chinrest and viewed 
the printed object next to a rendering of it on a color-calibrated dis¬ 
play. The color characterization of the display was performed from 
the chinrest position considering the viewing booth illumination. 
In this way, tristimulus values from objects placed on the display 
and illuminated by the viewing booth could be reproduced on the 
screen if viewed from the chinrest position. 
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0.269 0.322 0.372 0.409 


Fig. 13. Increasing number of layers results in increasing gamut volume. From left to right: 3, 6, 12 and 18 layers. The gamut volume is indicated below each 
image as a fraction of the volume of sRGB gamut shown transparent. 


Experimental Setup 



Viewing booth 

[LED Color Viewing Light M] 

• LED simulated CIED50 illuminant 

• Diffuse illumination perpendicular 
to the display's screen plane 


Chin rest 

Ensures: 

• ~60 ppd at the screen center 

• ~45® to the surface normal 

Display 

[Eizo ColorEdge cgSOlw] 

• 30" color calibrated screen 


Fig. 14. Physical setup for our visual experiment. 
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